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Abstract 

The linearized governing equations for an ideal fluid are presented for numerical analysis 
for the stability of free capillary surfaces in rectangular containers against unfavorable 
disturbances (accelerations, i.e. Rayleigh-Taylor instability). The equations are solved for 
the case of the right square cylinder. The results are expressed graphically in terms of a 
critical Bond number as a function of system contact angle. A critical wetting phenomena 
in the corners is shown to significantly alter the region of stability for such containers in 
contrast to simpler geometries such as the right circular cylinder or the infinite rectangular 
slot. Such computational results provide additional constraints for the design of fluids 
systems for space-based applications. 

INTRODUCTION 

Particularly since the inception of space flight a number of studies have been conducted 
to identify the stability limits of capillary surfaces to unfavorable disturbances (accelera- 
tions). The motivation for such investigations is generally to obtain design characteristics 
and performance limitations for in-space fluids management systems. For example, it is 
essential to understand the potentially destabilizing effects of a thruster firing on the liq- 
uid fuel in a partially-filled tank. In this paper a brief review of interfacial stability of 
the Rayleigh-Taylor-type will be provided which focuses on the restricted set of container 
geometries for which solutions are offered in the literature. In light of the growing need 
for design specific solutions for an ever increasing number of fluids systems applications, 
i.e. fluids experiments in space (Singh 1996), a new problem is outlined and solved for 
the stability of capillary surfaces in containers of rectangular cross-section. The results of 
this investigation, presented in terms of a critical Bond number as a function of contact 
angle and easily extendable to include container aspect ratio, may be readily added to the 
repertoire of the space systems designer. 

Review 

Surface tension forces dominate fluid interface behavior for low Bond number systems, 
B<1, where 
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Table 1: Correlation constants for eq. 3. 


where p is the density difference across the interface, a is the surface or interfacial tension, R 
is a characteristic dimension of the container, and g is the acceleration field strength. As B 
approaches 0(1), however, a critical balance is reached and, depending on the orientation 
of the acceleration field, further increases in g can cause destabilization of the interface 
and reorientation of the fluid to a perhaps undesirable location within the container. The 
precise value of B at which such a “reorientation” might occur (= Bo-) is an important 
design parameter for any capillarity-controlled fluid system and is particularly significant 
for fluids management processes in space. 

Critical Bond number analyses for confined geometries were performed indirectly as early 
as Duprez (1854) and Maxwell (1890) for the stability of pinned interfaces in circular and 
rectangular containers. These investigations are restricted to predominately flat surfaces 
originating out of an assumption of either a 90° contact angle condition or a pinned contact 
line. Treating the contact angle as a parameter and thus allowing significant curvature of 
the interface, as is most common in capillary systems, solutions were obtained by Concus 
(1968) for the right circular cylinder, Concus (1963) for the infinite slot, and Seebold et al. 
(1967) for the circular annulus. The experimental works of Masica et aL (1964) and Masica 
and Petrash (1965) concerning the cylinder and Labus (1969) concerning the annulus are 
also noteworthy. Solutions for spherical containers are presented by Reynolds and Satterlee 
(1966) and a number of solutions are reported for semi-bounded surfaces such as wall bound 
drops and bubbles by Reynolds and Satterlee (1966), pedant drops by Wente (1980), and 
liquid bridges by Coriell et al. (1977). Unbounded liquid layers are treated by Yiantsios 
and Higgins (1989) with pertinent references contained therein. 

Numerical correlations for Bo- as a function of the contact angle may be derived from 
the work of Concus (1963 and 1968) and Seebold et al. (1967), respectively. For the infinite 
rectangular slot one finds 

Bo- = 0.71 + 1.74 sin 0, (1) 

where B is defined using the slot half-width, and for the cylinder 

Bcr = 0.81 +2.59 sind, (2) 

where Bcr is defined using the cylinder radius. Larger values B > B^ lead to instability 
and breakup of the interface. Similar results for the annulus may be obtained. However, 
these depend heavily of the radius ratio Ri/Ro, where Ri and Ro are the inner and outer 
radii, respectively. In Fig. 1, correlations for B^ for free annular surfaces are plotted versus 
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FIGURE 1. Correlations of B „ with 6 for annular interfaces for a variety of radius ratios, 
Ri/Ro (numerical data of Seebold et al. (1967)). 

contact angle for a variety of radius ratios. The curves are determined using the form 

Bcr = Ci + Ci sin 0, (3) 

and the correlation constants C\ and C 2 are listed in Table 1. Note that Bcr is again 
defined on the outer radius, Ro, and that Ri/Ro = 0 recovers the correct form of eq. 3 for 
the circular cylinder, eq. 2. Observation of eqs. 1-3 and Fig. 1 reveals that B c r is nonzero 
and positive for all values of the contact angle 0. In addition, the solutions are symmetric 
about 0 = 90°. 

The commonality between circular cylinders, slots, and annular containers are the 
smooth continuous boundaries within which the fluids are confined. For the case of a 
container with an interior comer, the situation is altered significantly. As mathematically 
demonstrated by Concus and Finn (1969), when 0 < 90° - a (or 0 > 180° - a), here- 
after referred to as the Concus-Finn condition, a critical wetting condition is established 
resulting in complete wetting of the comer by the fluid: the fluid is pumped into and along 
the interior corners of the container by capillary forces, a is the corner half-angle. Such 
surfaces are unconditionally unstable to adverse accelerations. Thus, for fluid-container 
systems satisfying the Concus-Finn condition, Bcr = 0. Therefore, for problems such as 
the stability of a capillary surface in a rectangular container where a = 45°, nonzero values 
for Bo- may only be obtained for contact angles in the range 45° < 0 < 135°. A sketch of 
the different interfacial regimes is provided in Fig. 2 for a container of square cross-section. 
The condition of Fig. 2b (45° < 0 < 135°) is investigated here since the cases of Fig. 2a 
and 2c are unconditionally unstable for B > 0. 

ANALYTICAL SOLUTION 

Maxwell (1890) predicted the stability of an inverted capillary surface in a rectangular 
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(a) (b) (c) 

FIGURE 2. Sketch of wetting regimes in a square cross-sectioned container, liquid is shaded 
(note: container is bisected across diagonal): a. 0 < 45°, b. 45° < 6 < 135°, c. 6 > 135°. 


container of half-length a and half-width 6. Fig. 3 depicts the geometry under discussion, 
where the more dense fluid is below the interface and g acts positive in the positive in- 
direction. Maxwell’s solution is derived by minimizing the surface-plus-gravitational energy 
and assumes a pinned, predominately flat interface ( 6 « 90°). His result may be cast in 
terms of J3cr such that 



where 0 < b/a < 1. This solution approach may be extended to the case of perfect slip at 
the contact line where 0 is fixed at 90°. The result is 



As is commonly observed in practice, comparison of eqs. 4 and 5 shows that stability is 
significantly enhanced by the pinned condition. Note also from eq. 5 that for b/a — ► 0, 
Bo- — *• 7r 2 /4 which is equivalent to Concus’ (1963) solution for the infinite slot for 6 = 90° 
(see eq. 1), as well as to the solution for unbounded liquid layers where the disturbance 
wavelength A = 46, Yiantsios and Higgins (1989). No further analytical solutions are 
possible which allow appreciable variation in 6. 

NUMERICAL SOLUTION 

The numerical solution to the idealized equations of fluid motion are overviewed below 
in a like manner to Concus (1963), the dimensions of the problem being extended to analyze 
the surface h = h(x, y, t). Conservation of mass leads to the 3-dimensional Laplace equation 
for the velocity potential, <j>. The kinematic condition is applied at the free surface and the 
pressure jump condition across the interface due to capillary forces is then incorporated into 
Bernoulli’s law for a transient, inviscid, incompressible, and irrotational fluid. The resulting 


4 






FIGURE 3. Rectangular solution domain, a x 6. 

second order nonlinear partial differential equation is subject to the contact angle condition 
along the container walls. At this point the equations are linearized and normal modes are 
introduced for the velocity potential and for a small perturbation r?(x, y, t) to the leading 
order static interface shape H(x,y, B). The numerical solution to the resulting governing 
equation requires solution to the eigenvalue problem via the evaluation of the determinant 
of the solution matrix for the disturbance r). A negative (positive) determinant implies 
growth of the disturbance and thus instability (stability). Bo- is determined by the zero 
value. 

The mathematical solution detail is provided below in dimensionless form. Lengths are 
scaled such that x ~ a, j/ ~ 6, z ~ b and the surface height as measured from the x-y plane 
is h ~ b. Note that e = b/a and that 0 < e < 1. Time is scaled by (b/g) 1 / 2 , pressure by a/6, 
and the velocity potential by (gb 3 ) 1 ^ 2 . Subscript notation for differentiation is employed 
throughout. The equations are posed for the full domain (|x| < 1, |y| < 1). 

Generalized Equations 

Continuity of mass leads to 

V 2 0 = 0, (6) 

subject to <p n = 0 on the container walls, where n is the normal to the wall. The kinematic 
condition on the free surface is 


C?(frxflx "f~ (frzflz — (fry fit — 0 


( 7 ) 


on z — h(x, y, t). The pressure jump condition due to curvature of the free surface is given 
by 

-P = V • — — = Ch, (8) 


(l + |V/i| 2 ) 


1/2 


where £ is an operator on h given by 

_ e 2 h xx (l + h 2 ) +hy y (l+ e 2 hl) - 2e 2 h x h v h xy 
Cl> ~ (l +e s ftI + ^) 3/2 

The contact angle condition along the container walls is given by 

n • k = cos 6, 


( 9 ) 
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where k is the inward normal to the walls and n is the outward normal to the free surface 
given by 

k = Ofl, 0, 0) along x = ±1, 

k = (0, =Fl» 0) along y = ±1, 

and 

n = (l 4- t 2 h\ + hty 1 (-chx, -hy, 1) , 

respectively. Thus, the contact angle conditions at the boundary of the surface are 

± — rrr = CO6 0 along X = ±1 (10) 

(l + e^f + fc*)' 

and 

±- — —rjo = cos 6 along y = ± 1. (11) 

(l+eW x + hl) 

Incorporating eq. 8 into Bernoulli’s equation for a transient, ideal fluid yields 

B + - (t 2 <t>x + + hj + Ch — C, (12) 

which is applicable on z = h(x, y, t). C in the above equation is a constant, and in the most 
general sense C = C(t) and is determined by the volume of fluid present in the container, 
which is here assumed steady in time. 

Linearized Governing Equations 

Introducing the perturbation 


h = H(x,y)+r)(x,y,t), 

normal modes are selected for </> and 77 such that 

(13) 

<f> = <p'{x, y, z ) cos(u><t), 

(14) 

77 = rf(x,y) sin(a;jt). 

(15) 


Substituting eqs. 13-15 into eq. 12, neglecting nonlinear terms, and noting Ui = 0 for 
neutral stability, yields the simplified Bernoulli equation 


B(H + t)) + £(H + ti) = C, (16) 

where the prime notation for 77 has been dropped for clarity. Assuming rj/H 1, the 
zeroeth order solution for the interface shape H may be determined from eq. 16 to be 

BH + CH = C, (17) 
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where CH is the operation of eq. 9 on H. Eq. 17 is subject to 
. eH x 


and 


(l+e 2 //| + ff y 2 ) 1/2 

Hz 

.1/2 


= cos 6 along x = ±1 


= cos# along y = ±1. 


(l + e 2 H 2 + H 2 ) 

The first order solution for the perturbation rj is given by 

Brj + Crj = 0, 

subject to 


(18) 


erj x (l + Hy'j — eT] y H x Hy = 0 along 


x = ±1 


and 


c n = 


% (l + e 2 H f) - t 2 T) x H x H y = 0 along y = ±1. 

The operation Crj is defined by 

7 U 2 Tlxx (i + Hi) + rjyy (l + e 2 Hl) - 26 2 7 lxy H x H y 

(l+e 2 H 2 + H 2 )' 1 V ' V J 

“f~2£ {flyHxxHy ~f _ TjxHxHyy TJyHxHxy — 7^x*Hy-Hzy)j 


(19) 


and is determined by expanding C(H + rj) in powers of rj retaining only terms of 0{rj). In 
the limit £ 2 « 1, 

Vyy ^TjyHyHyy 

* 3/2 / . _\ 5 / 2 » 


Crj = 


(i (•+«?)' 

subject to the leading order boundary condition r\y = 0 on y = ±1, which recovers the 
governing system of Concus (1963) for the infinite slot. It is important to note that for this 
limiting case the first boundary condition to eq. 18 is O(e) and is ignored. The infinite slot 
formulation is thus accurate to 0(e) for a finite slot provided 6 90°. 


Numerical Solution Detail 


In the numerical solution procedure eq. 17 is discretized based on a fourth-order central- 
differencing scheme. For the calculation of the static shape of the free surface H(x,y;B), 
the Newton iteration method with successive under-relaxation is used to address the non- 
linearity of the governing equation. In determining the eigenvalue, B = Bo-, the same 
discretization as that used for solving the static interface shape is used and can be ex- 
pressed in the following form 

[A + BI] rj = 0. (20) 
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FIGURE 4. Surfaces profiles across container diagonal of a square container for a variety 
of contact angles and Bond numbers. 


The determinant of the coefficient matrix in eq. 20 is zero at B = B a •. Since the coefficient 
matrix [A] is dependent on the static shape and thus B, the overall solution procedure 
involves an iteration between eq. 20 via eq. 18 and the calculation of the static shape, 
eq. 17. A bisection method is used to determine B „ in the iteration process. It should be 
noted that the determination of the critical Bond number has to be based on the solution 
from the full domain (i.e. all four solid boundaries included). This is due to the fact that 
an asymmetric disturbance leads to the fundamental subharmonic mode instability. 

As 0 decreases towards 45° (or increases towards 135°), the nonlinearity of the static eq. 
17 increases dramatically, requiring a smaller relaxation factor for the Newton iterations. 
Hence, the run time per case increases as the contact angle deviates from 90°. For a 60 x 60 
grid system, the typical run time per case on an SGI Indigo II with a single-processor is one 
hour for contact angles in the vicinity of 90°. However, the run time becomes significantly 
longer, reaching 24 hours, when the contact angle is close to 45°. Expectedly, solutions for 
Bcr are found to be symmetric about 6 = 90°. 

RESULTS 

In Fig. 4, surface profiles across the diagonal of a square cross-section container (e = 1) 
are compared for a variety of Bond numbers for contact angles 86° and 60°. The coordinates 
for the figure are normalized by the diagonal of the container and are presented to scale. The 
base-state interface shape H(x,y) for the case e = 1, 6 = 60° is presented 3-dimensionally 
in Fig. 5 for B = 3.0. Slight inflections of the interface near the corners are observed 
which can also be discerned in Fig. 4, 6 = 60°, B „ = 3.0. These appear to diminish with 
decreasing 6. 
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FIGURE 5. Surface shape H(x,y;Ba.) for 9 = 60° in a square cross-sectioned container, 
e = 1. 


The numerical results for B^ are presented in Fig. 6 as a function of 6 for the aspect 
ratio 6 = 1. The numerical results of Concus (1963) for the infinite slot and the right circular 
cylinder are also provided via eqs. 1 and 2. The region of stability denoted by the area 
below the curves is obviously altered for the rectangular section of this study when compared 
to the infinite slot. This is attributable to the restricted range of contact angles allowing 
for stable interfaces which cover the solution domain. Further “preliminary” calculations 
show that the rapid decrease of Ba- towards zero for 9 48° reveals the sensitivity of the 

surface to the critical contact angle condition which is satisfied for 9 < 45°. This trend is 
indicated on Fig. 6 by the dashed line extrapolations from the numerical solutions obtained 
for 60° < 9 120°. It is useful to note that for 9 £ 48°, B takes normally anticipated 

values, 0(1). The effect of contact angle hysteresis on such stability results is likely to 
delay the instability while equilibrium conditions are established at the contact line. If 
the disturbance has temporal periodicity, the effect of hysteresis could be to significantly 
increase the stability of the interface, particularly for contact angles near 45°. However, as 
found in recent space experiments by Concus et al. (1997), extended periods of thermal 
and mechanical disturbances in the presence of a steady background acceleration such that 
B £ Bcr will ultimately bring about the predicted instability. 

It is important to note that the rectangular geometry of this investigation is funda- 
mentally different from the infinite rectangular slot as seen by the limiting case of e — ► 0. 
One might expect the solution to agree in this limit, however, the presence of the corners 
dramatically alters the base state surface profile for 9 \ 45°, or 9 f 135°. 

It is also of interest to note that an inflection point appears in the static interface 
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FIGURE 6. Bcr versus 0 for e = 1; ■ numerical results, — eqs. 1 and 2 for the infinite 
rectangular slot and right circular cylinder, and y and □ are eqs. 4 and 5 for the case 
6 = 1,0 = 90°, respectively. Dashed line is extrapolation of numerical results. 


shape as a result of the additional space dimension in the solution of the problem in x 
and y. Such inflections are not found in the cases of the slot and the cylinder which are 
1-dimensional problems where the onset of any inflection of the surface actually signals the 
onset of instability, Concus (1964). 

The fact that Bcr determined numerically agrees with eq. 5 (y on Fig. 6) in the limit 
6 -» 90° is expected and serves in part as a verification of the numerical results. The code 
was also “checked” using 90° contact angle conditions on opposing faces of the containers 
while varying 0 on the other opposing faces. The results from these runs recovers the infinite 
slot results of Concus (1963) for all values of the contact angle. For the full numerical 
problem, however, for contact angles lower than « 60°, the numerical approach employed 
experiences convergence difficulties. What is observed is that for 0 < 60° (0 > 120°) 
the L <2 norm (= ||if n+1 — 7f n ||) of the base state eq. 17 decreases to a minimum and 
then proceeds to increase after excessive iterations without achieving machine zero. It is 
taken for granted that the steepening of the interface in the corner as the contact angle 
decreases, typified by the angle of the interface in the corner measured along the corner 
bisector 0 C = cos -1 (\/2cos0), leads to increased grid resolution uncertainties, increased 
nonlinearities in the governing equations, and thus increased run time and roundoff errors. 
But these effects should not be appreciable for say 0 55°, where 0 C = 41°. Since a slight 

variation in B^- was also detected when 0 £ 60° for different values of the relaxation factor, 
it is our suspicion that the Newton iteration method with under-relaxation may not be the 
most robust technique for the system of equations solved herein. As a recourse, a time- 
marching scheme is presently being developed to further investigate the numerical issues. 1 
Regardless, the qualitative nature of the computational results will remain unchanged as 
those presented in Fig. 6. With the numerical issues resolved, calculations for Bo- will be 

1 The dashed line extrapolations of Fig. 6 are actually numerical solution curves deemed “pre- 
liminary” in that these experience convergence difficulties. 
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completed for a range of aspect ratios, e. 

CONCLUSION 

The governing equations and boundary conditions for the determination of the dynamic 
stability of capillary surfaces in rectangular containers are extended to 3-dimensions and 
presented for numerical solution. Calculations for a square cro6s-sectioned container are 
performed which serve as a model for containers possessing interior corners. Such container- 
types are commonly employed in fluids management systems in space. The results reveal 
that, though stability is comparable to the circular cylinder for large contact angles near 
90°, the range of contact angles yielding positive values for the critical Bond number is 
significantly reduced due to a corner wetting phenomena governed by the Concus-Finn 
condition. Bcr is determined to be 0(1) for 48° < 0 < 132°, but diminishes rapidly to zero 
as 6 approaches 45° from above, or 135° from below. Computations addressing the effects 
of container aspect ratio e are currently underway. 
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